void Gauss_Elimination_Row(std::vector<std::vector<double> > &_A, 
                           std::vector<double> &_b);
void Gauss_Elimination_Row_Main(std::vector<std::vector<double> > &_A, 
                                std::vector<int> &_u);



void Gauss_Elimination_Row(std::vector<std::vector<double> > &_A, std::vector<double> &_b)
{
    int n = _A.size();
    
    std::vector<std::vector<double> > L(n, std::vector<double> (n, 0));
    std::vector<std::vector<double> > U(n, std::vector<double> (n, 0));
    std::vector<int> u(n, 0);
    
    Gauss_Elimination_Row_Main(_A, u);
    Get_L(_A, L);
    Get_U(_A, U);
    int i = 0;
    double tmp;
    while(i <= n - 2)
    {
        if(i != u[i])
        {
            tmp = _b[i];
            _b[i] = _b[u[i]];
            _b[u[i]] = tmp;
        }
        ++i;
    }
    Calculate_LU(L, U, _b);
    


    for(i = 0; i < n; ++i)       //change -0 to 0
    {
        if(_b[i] == -0)
        {
            _b[i] = 0;
        }
    }
}

void Gauss_Elimination_Row_Main(std::vector<std::vector<double> > &_A, std::vector<int> &_u)
{
    int i, j ,k;
    int n = _A.size();
    int p;
    for(k = 0; k <= n - 2; ++k)
    {
        //find p
        p = k;
        for(i = k + 1; i <= n - 1; ++i)
        {
            if(fabs(_A[i][k]) > fabs(_A[p][k]))
            {
                p = i;
            }
        }

        if(p != k)
        {
            Exchange_Line(_A, p, k, 0, n - 1);
        }
        _u[k] = p;
        if(_A[k][k] != 0)
        {
            for(i = k + 1; i <= n - 1; ++i)
            {
                _A[i][k] = _A[i][k] / _A[k][k];
                for(j = k + 1; j <= n - 1; ++j)
                {
                    _A[i][j] = _A[i][j] - _A[i][k] * _A[k][j];
                }
            }
        }
        else
        {
            break;
        }
    }
}

void Gauss_Elimination_All_Main(std::vector<std::vector<double> > &_A, std::vector<int> &_u, std::vector<int> &_v)
{
    int i, j ,k;
    int n = _A.size();
    int p, q;
    for(k = 0; k <= n - 2; ++k)
    {
        //find p
        p = k;
        for(i = k + 1; i <= n - 1; ++i)
        {
            if(fabs(_A[i][k]) > fabs(_A[p][k]))
            {
                p = i;
            }
        }

        //find q
        q = k;
        for(i = k + 1; i <= n - 1; ++i)
        {
            if(fabs(_A[k][i]) > fabs(_A[k][q]))
            {
                q = i;
            }
        }
        

        if(p != k)
        {
            Exchange_Line(_A, p, k, 0, n - 1);
        }
        _u[k] = p;

        if(q != k)
        {
            Exchange_Row(_A, p, k, 0, n - 1);
        }
        _u[k] = p;


        if(_A[k][k] != 0)
        {
            for(i = k + 1; i <= n - 1; ++i)
            {
                _A[i][k] = _A[i][k] / _A[k][k];
                for(j = k + 1; j <= n - 1; ++j)
                {
                    _A[i][j] = _A[i][j] - _A[i][k] * _A[k][j];
                }
            }
        }
        else
        {
            break;
        }
    }
}

void Gauss_Elimination_All(std::vector<std::vector<double> > &_A, std::vector<double> &_b)
{
    int n = _A.size();
    
    std::vector<std::vector<double> > L(n, std::vector<double> (n, 0));
    std::vector<std::vector<double> > U(n, std::vector<double> (n, 0));
    std::vector<int> u(n, 0);
    std::vector<int> v(n, 0);
    
    Gauss_Elimination_All_Main(_A, u, v);
    Get_L(_A, L);
    Get_U(_A, U);
    int i = 0;
    double tmp;
    while(i <= n - 2)
    {
        if(i != u[i])
        {
            tmp = _b[i];
            _b[i] = _b[u[i]];
            _b[u[i]] = tmp;
        }
        ++i;
    }
    Calculate_LU(L, U, _b);
    for(i = n - 2; i >= 0; --i)
    {
        if(i != v[i])
        {
            tmp = _b[i];
            _b[i] = _b[v[i]];
            _b[v[i]] = tmp;
        }
        ++i;
    }


    for(i = 0; i < n; ++i)       //change -0 to 0
    {
        if(_b[i] == -0)
        {
            _b[i] = 0;
        }
    }
}
